Zonal Flow as Pattern Formation: Merging Jets and the Ultimate Jet Length Scale 



Jeffrey B. Parker* and John A. Kromnies' 
Princeton University, Princeton Plasma Physics Laboratory, Princeton, New Jersey 08543, USA 

(Dated: January 23, 2013) 

Zonal flows are well known to arise spontaneously out of turbulence. It is shown that for statisti- 
cally averaged equations of quasigeostrophic turbulence on a beta plane, zonal flows and inhomoge- 
neous turbulence fit into the framework of pattern formation. There are many implications. First, 
the zonal flow wavelength is not unique. Indeed, in an idealized, infinite system, any wavelength 
within a certain continuous band corresponds to a solution. Second, of these wavelengths, only those 
within a smaller subband are linearly stable. Unstable wavelengths must evolve to reach a stable 
wavelength; this process manifests as merging jets. 



Introduction Zonal flows — azimuthally symmetric, 
generally banded, shear flows — are spontaneously gen- 
erated from turbulence and have been reported in atmo- 
spheric [1], oceanic [2], and laboratory plasma [3] con- 
texts. Recently, they have also been observed in astro- 
physical simulations [4]. In magnetically confined plas- 
mas, zonal flows are thought to play a crucial role in 
regulation of turbulence and turbulent transport [5]. 

Zonal flows remain incompletely understood, even re- 
garding the basic question of the jet width. Various au- 
thors have attempted to relate the jet width or spacing 
to length scales that emerge from the vorticity equation 
by heuristically balancing the magnitudes of the Rossby 
wave term and the nonlinear advection term. Those 
scales include the Rhines scale (J7//3) 1 / 2 , where U is the 
rms velocity and (3 is (in planetary contexts) the north- 
ward gradient of the Coriolis parameter [6], and a re- 
lated transitional scale obtained by using the energy in- 
put £ rather than U [7]. It is often asserted that an 
inverse energy cascade is arrested at one or another of 
those scales by the generation of Rossby waves, which 
then preferentially sends energy into zonally elongated 
structures. However, others have concluded that the in- 
verse cascade cannot be stopped by a f3 effect but only 
by friction [1, 8]. Furthermore, it has been suggested 
that the primary mechanism maintaining the jets is in- 
stead a spectrally nonlocal interaction with turbulence at 
small scales [9] . Regardless of the physical mechanism in- 
volved, heuristic balances are not systematically derived 
or quantitatively predictive. 

A related topic is the merging of jets. Coalescence 
of two or more jets is ubiquitous in numerical simula- 
tions [9-11]. The merging process occurs during the ini- 
tial transient period before a statistically steady state is 
reached. It is clear that the merging is part of a dy- 
namical process through which the zonal flow reaches its 
preferred length scale, but there has been little theoreti- 
cal understanding of the merging phenomenon. Insights 
into jet structure are valuable for explaining various fea- 
tures of weather and for untangling a host of nonlinear 
processes in plasmas, including details of transitions be- 
tween modes of low and high confinement. 

Our present work addresses these questions in the con- 



text of stochastically forced quasigeostrophic (QG) flow 
on a P plane, a model for fluid turbulence in a rotating 
system [7]. It is known that the QG equation for po- 
tential vorticity is nearly mathematically identical to the 
generalized Hasegawa-Mima equation [12, 13] for electro- 
static potential, a model for magnetized plasma turbu- 
lence in the presence of a background density gradient. 
Importantly, numerical simulations of both models can 
display emergence of steady zonal flows. We concentrate 
on the QG model because it has been studied far more 
extensively in the literature. However, our methodology 
is equally applicable to the Hasegawa-Mima equation, as 
we will describe in a future, more detailed article. 

We study a statistical average of the QG equation. Sta- 
tistical approaches enable one to gain physical insight by 
averaging away the details of the turbulent fluctuations 
and working with smoothly varying quantities. Some- 
times, statistical turbulence theories strive for quantita- 
tive accuracy, which requires rather complicated methods 
[14]. In contrast, our investigation is at a more basic level 
and concerns the fundamental nature of zonal flows inter- 
acting sclf-consistently with inhomogeneous turbulence. 

We discover that the statistically averaged system pos- 
sesses the mathematical structure of pattern formation. 
Pattern formation is the study of systems out of thermal 
equilibrium that spontaneously generate spatial patterns 
[15-20]. Application of its theory to the averaged system 
leads to insights that are difficult to discern by studying 
one realization of the turbulence. 

Two important results follow that are general proper- 
ties of pattern formation systems. First, the zonal flow 
wavelength is not unique. Indeed, in an idealized, infi- 
nite system, any wavelength within a certain continuous 
band corresponds to a steady-state solution. Second, of 
these wavelengths, only those within a smaller subband 
are linearly stable. Unstable wavelengths must evolve to 
reach a stable wavelength. For short (long) wavelength 
unstable jets, this process manifests as merging (branch- 
ing) jets. 

Quasigeostrophic equation and CE2 model Our basic 
model is two-dimensional quasigeostrophic turbulence on 
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a (3 plane in the limit of infinite deformation radius, 

d t w + v • Vw + I3d x <j) = £-nw- u(-l) h V 2h w, (1) 

where w is the relative vorticity, <f> is the streamfunction 
such that w — V 2 0, v = z x V0 is the horizontal fluid 
velocity, i is in the vertical direction, /i is the constant 
friction, v is the viscosity with hyperviscosity factor h, 
and £ is small-scale white-noise forcing. The zonal flow 
behavior in numerical simulations of Eq. (1) is shown in 
Fig. 1(a). During the transient period, merging jets are 
observed, while in the late time, a statistically steady 
state is reached with stable unwavering jets. 

We restrict ourselves to the quasilinear (QL) approxi- 
mation of this system. To obtain the QL equations, we 
perform an eddy-mean decomposition, given by decom- 
posing all fields into a zonal mean and a deviation from 
the zonal mean, then neglect the eddy-eddy nonlinear- 
ities within the eddy equation [11]. The QL model ex- 
hibits the same basic zonal jet features as the full model, 
namely merging jets and the formation of stable jets. 

We consider a statistical average of the QL system. In 
the presence of steady zonal flows, a statistical homogene- 
ity assumption is clearly invalid. Therefore, we allow the 
turbulence to be inhomogeneous in the direction of zonal 
flow variation (y). The averaged equations, referred to 
as the second-order cumulant expansion (CE2), are [11] 



(a) 



200 



dtW 



{U 1 -U 2 )d x W -{U'{ -U 2 r ) 
- [2/? - {U'-l + U'i)\dyd v d x C 

= F — 2fiW - 2vD h W, (2a) 
d t U + dyd y d x C(0,0,y,t) = -fOJ - u{-\) h dfU, (2b) 

where x and y represent two-point separations, y rep- 
resents the two-point average position (if the turbu- 
lence were homogeneous, there would be no y depen- 
dence), W(x,y | y, t) and C(x,y | y,t) are the one- 
time, two-space-point correlation functions of vorticity 
and streamfunction, U (y, t) is the zonal flow velocity, 
C7i = U(y + y/2,t), U 2 = U(y-y/2,t), F(x,y) is chosen 
to be isotropic, homogeneous ring forcing as in [11], and 
Dh is the hyperviscosity operator. The relation between 
W and C can be found in [11]. 

Given the assumption that the stochastic forcing £ is 
white (delta-correlated) noise, the only further assump- 
tions necessary for CE2 to be an exact description of the 
QL model are statistical homogeneity and ergodicity in 
the zonal (x) direction. This is because the QL model 
neglects the nonlinear eddy-eddy term that would give 
rise to a closure problem. Alternatively, CE2 can be re- 
garded as a drastically truncated statistical closure of the 
full QG model [21-24]. However, for present purposes we 
prefer the former interpretation. 

The CE2 equations exhibit several important symme- 



»» 100 


(b) 40 
20 

s» 
-20 
-40 









50 



100 

t 



150 



200 



FIG. 1. (a) Merging jets during the transient regime of the 
QG equation (1) (zonal-mean velocity is shown), (b) Merging 
behavior observed in the amplitude equation (4) [KeA(y, t) is 
shown] . 



tries of translation and reflection, given by 

y^y + Sy, (3a) 
x,y -> -x, -y, (3b) 

y,y^-y,-y, (3c) 

x,y-*-x,-y. (3d) 

CE2 has been studied extensively numerically [21-24]. 
Those simulations also exhibit merging jets [22] . 

Zonostrophic Instability and the Amplitude Equation 
For Eq. (2) there always exists a homogeneous equilib- 
rium, which arises from a simple balance between forcing 
and dissipation: W = (2/i + 2vD h )- 1 F, U = 0. This 
equilibrium is stable in a certain regime of parameters. 
As a control parameter such as the friction fi is varied, 
this homogeneous state becomes zonostrophically unsta- 
ble [11, 22]. Physically, zonostrophic instability occurs 
when dissipation is overcome by the mutually reinforcing 
processes of eddy tilting by zonal flows and production 
of Reynolds stress forces by tilted eddies. The instabil- 
ity eigenmode consists of perturbations spatially periodic 
in y with zero real frequency [11], so that zonostrophic 
instability arises as a Type I s instability [15] of homoge- 
neous turbulence. 

Just beyond the instability threshold, a bifurcation 
analysis yields a perturbative amplitude equation for 
the bifurcating mode. This amplitude equation is con- 
strained by the translation and reflection symmetries (3) 
to take a universal form [15]. The amplitude equation, 
sometimes referred to as the real Ginzburg-Landau equa- 
tion, is 



d t A{y,t)=A + dU-\A\ 2 A, 



(4) 



where all coefficients have been rescaled to unity. Here, A 
is the complex, spatially varying amplitude of the eigen- 
vector that is neutrally stable at the bifurcation point. 
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The minus sign in the last term corresponds to a super- 
critical bifurcation, which has been demonstrated numer- 
ically to occur [22]. 

The amplitude equation (4) is well understood [15-17], 
and much of its qualitative behavior is seen generically 
in pattern formation systems. First, a steady-state so- 
lution exists for any wave number within the continuous 
band — 1 < k < 1 (to see this, observe that A = ae tky 
with |a| 2 = 1 — k 2 is a solution). Second, only solu- 
tions with k 2 < 1/3 are stable [16]. This is demonstrated 
in Fig. 1(b), where an unstable solution that has been 
slightly perturbed undergoes merging behavior until a 
stable wave number is reached. 

These qualitative behaviors are also exhibited by the 
CE2 system. In the following sections we calculate the 
equilibria and stability of CE2. 

Calculation of Ideal States We proceed to find the 
steady-state solutions of Eq. (2). In the context of an 
infinite domain with no boundaries, these solutions are 
referred to as ideal states. Let q denote the basic zonal 
flow wave number of an ideal state. For a given q, we 
solve the time-independent form of Eq. (2) directly. This 
approach is distinct from time integration of Eq. (2) to 
a steady state, which is done in [21-24] within a finite 
spatial domain. Our procedure has two advantages, both 
related to the fact that ideal states exist for any q within 
a continuous band. First, we can specify precisely the q 
of the desired solution. Second, we can solve directly for 
all solutions, including unstable ones, rather than find 
only those which develop from time evolution. 

An ideal state is represented as a Fourier-Galerkin se- 
ries with coefficients to be determined [16, 19, 20]. We 
expand as follows: 



U(y) = £ U p e**, 
p=-p 

MNP 

W(x,y\y) = E E W mnp e ma *e mb y e 

m=-M n=—N p= — P 



(5a) 



ipqy 



(5b) 



While the periodicity in y is desired, the correlation func- 
tion should decay in x and y; periodicity in x and y is a 
consequence of using the convenient Fourier basis. Thus, 
a and 6, unlike q, are numerical parameters. They rep- 
resent the spectral resolution of the correlation function 
and should be small enough to obtain an accurate solu- 
tion. 

The CE2 symmetries allow us to seek a solution where 
U(y) = U(-y) and W(x,y | y) = W(-x,-y | y) = 
W(x,~y | — y) = W(~x,y \ —y). These constraints, 
along with reality conditions, force U p to be real, U p — 



U- p , and W mnp = W*_ 



= W* 



= W* 



mnp Tr —m,Ti,p " m,—n,p ' ' m,n } —p' 

We obtain a system of nonlinear algebraic equations for 
the coefficients U p , W mnp by substituting the Galerkin 
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FIG. 2. Zonal flow amplitude Ui, U2 as a function of ideal 
state wave number q at (a) fx = 0.21 and (b) fi = 0.19. In the 
unshaded region, ideal states are stable. The vertical lines 
correspond to various instabilities which separate the regions 
(see Fig. 3). 



series into Eq. (2) and projecting onto the basis func- 
tions. To demonstrate the projection for Eq. (2a), let 
4>mnp = e lmax e mby e l P q y. We project Eq. (2a) onto (j> rst 
by operating with 



2 7 r2 7 r2 7 rV 1 j , T 7 " _ 
7 / dx dy dy<j> rst . 

a q J J-tt/o, J-ir/b J-iv/q 



(6) 



For instance, the term (U\ — U-2)d x W projects to 
Up' W mnp , where repeated indices are summed 



J (1) , 

rstp'mnp 

0Ver ' ^llp'rnnp = ima{ai - (T2)S m ,r$p>+ P -t, 01,2 = 

sinc(ai.27r/&), and ai^ = nb — sb ± p'q/2. The other 
terms of Eq. (2a), as well as Eq. (2b), are handled sim- 
ilarly. In total, we generate as many equations as there 
are coefficients. 

The system of nonlinear algebraic equations is solved 
with a Newton's method [25]. Figure 2 shows the zonal 
flow amplitude coefficients U p as functions of q at /i = 
0.21 and ijl = 0.19. Near the instability threshold, ideal 
states exist at all q for which the homogeneous equilib- 
rium is zonostrophically unstable [between the two lines 
labeled N in Fig. 2(a)]. Farther from threshold, there 
is a region of q where the ideal state solution disappears 
[between the lines N and D in Fig. 2(b); see also Fig. 
3]. The values of the other parameters used are (3 = 1, 
v = 10~ 3 , and h = 4. The forcing F acts equally on wave 
vectors k with 7/8 < |k| < 9/8. The forcing provides a 
total energy input of e = 1. 

Stability of the Ideal States To investigate stability of 
the ideal states, we consider perturbations SW(x, y \y,t) 
and SU(y,t) about an equilibrium W,U. The linearized 
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CE2 equations are 

d t SW = - [SUi ~ 5U 2 )d x W - (U! - U 2 )d x SW 
+ (6U[' - 6U%)(V 2 + \£%)d x C 
+ (U'{ + Etf) (V 2 + l4)d x 6C 
- {SU'{ + 5U' 2 r )dyd y d x C - {U'( + U^)dyd v dJC 
+ 2/3dyd y d x 6C - 2(/x + vD h )5W, (7a) 

d t su = -[n + v{-i) h df] su - dydydjc (0,0 \%t). 

(7b) 

Since the underlying equilibrium is periodic in y, the per- 
turbations can be expanded as a Bloch state [16, 19]: 

SW(x,y\ y,t) = e at e lQ y ^ SW mnp e max e mby e lpqV , 

mnp 

(8a) 

6U(y, t) = e at e lQ y ^ 5U p e ip ™ \ (8b) 
p 

where Q is the Bloch wave number and can be taken to 
lie within the first Brillouin zone —q/2 < Q < q/2. We 
do not use a Q x or Q y because as previously mentioned 
the periodicity in x and y is artificial. Equation (7) is 
projected onto the basis functions in the same way as in 
the ideal state calculation. This projection results in a 
linear system at each Q for the coefficients SW mnp and 
5U p ] this determines an eigenvalue problem for a. The 
equilibrium is unstable if for any Q there are any eigen- 
values with Re cr > 0. 

The stable ideal states exist within the closed region, 
known as the stability balloon, shown in Fig. 3. We adopt 
the friction /i as the control parameter, but plot using 
—\x because stability balloons are conventionally depicted 
with destabilization parameters. The stability balloon is 
bounded by curves representing marginal stability due to 
various instabilities. These curves are determined near 
threshold (/i = fi c = 0.24) by the Eckhaus instability, 
a long-wavelength universal instability, and farther from 
threshold by various short-wavelength instabilities. De- 
tails of these instabilities will be reported in a subsequent 
article. 

For 0.07 < (j, < fi c , the stability balloon is consistent 
with the QL zonal flow wave numbers obtained in [11]. 
On the other hand, for /i < 0.07 there is some discrep- 
ancy between CE2 and the QL system. CE2 predicts no 
steady states, and indeed, numerical simulations of the 
CE2 equations display translating zonal jets. In contrast, 
the QL system exhibits seemingly steady zonal jets for 
even for these small values of /i [11]. This discrepancy 
warrants further investigation. 

Discussion Numerical simulations typically occur 
within a finite domain. When periodic boundary con- 
ditions are used, our infinite-domain results are modified 
merely by the discretization of wave numbers. This af- 
fects not only the possible equilibria, but also any per- 
turbations and hence the stability boundaries too. 
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FIG. 3. Stability balloon for the CE2 equations. Above the 
neutral curve (N), the homogeneous turbulent state is zonos- 
trophically unstable. Ideal states are stable within the closed 
region. The solution branch of stationary ideal states van- 
ishes to the left of D. Near threshold the stability balloon is 
determined by the Eckhaus instability (E). 

For a time-evolving system, the exact q that is ulti- 
mately chosen within the stability balloon results from a 
dynamical process and is not addressed in a systematic 
way by the present study. 

While the CE2 equations exhibit spontaneously gener- 
ated zonal flows, it is true that they neglect many phys- 
ical effects. An important piece of physics missing from 
the CE2 equations is the nonlinear eddy self-interaction, 
which clearly cannot be ignored in general. Furthermore, 
the CE2 equations involve one-time correlation functions 
rather than the more general two-time functions. The 
lack of time-history information means that most of the 
effects of wave propagation arc discarded [26]. At least 
one particular instance of the qualitative failure of CE2 
has been noted [27]. 

Yet, the basic mathematical structure of the theory 
presented here arises only from symmetry arguments and 
general properties of the zonostrophic instability. If one 
were to include the important physics neglected in CE2, 
those general symmetries and properties should remain 
intact. Therefore, we expect our qualitative conclusions 
to likewise remain valid. 

In summary, by analyzing a second-order statistical 
model of an ensemble of interacting zonal flows and tur- 
bulence, we have shown that zonal flows constitute pat- 
tern formation amid a turbulent bath. We calculated the 
stability balloon of steady zonal jets and explained the 
merging of jets as a means of attaining a stable wave num- 
ber. In general, the use of statistically averaged equa- 
tions and the pattern formation methodology provides 
a path forward for further systematic investigations of 
zonal flows and their interactions with turbulence. 

We acknowledge useful discussions with Brian Farrell, 
Henry Greenside, Petros Ioannou, and Brad Marston. 
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